echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 9.10. (Finite-Difference Method).
% Section	9.9,	Finite-Difference Method, Page 496
echo on; clc; format long; hold off; clear

% This program implements the finite difference
% method for solving the boundary value problem

%       x`` = p(t)x`(t) + q(t)x(t) + r(t)    
% with  x(a) = alpha, x(b) = beta

%	Store p(t),q(t), r(t)	in the M-files p.m  q.m  r.m 
% function z = p(t)
% z = 2*t/(1 + t^2);

% function z = q(t)
% z = - 2/(1 + t^2);

% function z = r(t)
% z = 1;

delete p.m
diary p.m; disp('function z = p(t)');...
           disp('z = 2*t/(1 + t^2);');...
diary off;

delete q.m
diary q.m; disp('function z = q(t)');...
           disp('z = - 2/(1 + t^2);');...
diary off;

delete r.m
diary r.m; disp('function z = r(t)');...
           disp('z = 1;');...
diary off;

% Remark. p.m  q.m  r.m findiff.m trisys.m are used for Algorithm 9.10
p(0); q(0); r(0);
pause	% Press any key to continue.

clc;
%	Place the endpoints of [a,b] in  a  and  b.

%	Place the initial  value x(a) in  alpha.

%	Place the terminal value x(b) in  beta.

%	Place the number of subintervals in  n.

a  = 0;
b  = 4;
alpha =  1.25;
beta  = -0.95;
n  = 20;

[T,X] = findiff('p','q','r',a,b,alpha,beta,n);
points = [T;X];

pause	% Press any key to see the solution points.

clc; clg;
c = min(X)-0.2;
d = max(X)+0.2;
axis([a b c d]);...
plot(T,X,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
if n<=30,
  plot(T,X,'or');
end;...
xlabel('t');...
ylabel('x');...
Mx1 = 'The finite difference solution to x`` = f(t,x,x`).';...
title(Mx1);...
grid;...
axis;...
hold off;...
shg; pause	% Press any key to continue.

Mx2 = '     t(k)               x(k)'; clc;
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
diary off,echo on

pause	% Press any key to perform extrapolation.

%	Place the endpoints of [a,b] in  a  and  b.

%	Place the initial  value x(a) in  alpha.

%	Place the terminal value x(b) in  beta.

%	Place the number of subintervals in  n.

a  = 0;
b  = 4;
alpha =  1.25;
beta  = -0.95;
n  = 20;

[T,X1] = findiff('p','q','r',a,b,alpha,beta,n);
[T,X2] = findiff('p','q','r',a,b,alpha,beta,2*n);
X2 = X2(1:2:length(X2));
[T,X3] = findiff('p','q','r',a,b,alpha,beta,4*n);
X3 = X3(1:4:length(X3));
T = T(1:4:length(T));
Z1 = (4*X2-X1)/3;
Z2 = (4*X3-X2)/3;
X = (16*Z2-Z1)/15;
points = [T;Z1;Z2;X];

pause	% Press any key to see the solution points.

clc; clg;
c = min(X)-0.2;
d = max(X)+0.2;
axis([a b c d]);...
plot(T,X,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
if n<=30,
  plot(T,X,'or');
end;...
xlabel('t');...
ylabel('x');...
Mx1 = 'The finite difference solution to x`` = f(t,x,x`).';...
title(Mx1);...
grid;...
axis;...
hold off;...
shg; pause	% Press any key to continue.

format long
Mx2 = 'Extrapolated solutions z1(k), z2(k) and x(k).';
Mx3 = '     t(k)               z1(k)              z2(k)              x(k)';
clc,echo off,diary output,...
disp(''),
disp(''),disp(Mx1),disp(''),disp(Mx2),...
disp(''),disp(Mx3),disp(points'),...
diary off,echo on
